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TECHNICAL NOTE 4273 


ON PAIRS OF SOLUTIONS OF A CLASS OF INTERNAL VISCOUS 
FLOW PROBLEMS WITH BODY FORCES 
By Simon Ostrach and Lynn U. Albers 


SUMMARY 

In previous analyses of fully developed combined forced- and 
natural- convection flows, a few examples were presented of two distinct 
states of flow and beat transfer vbieh were obtained for a given set of 
conditions If the frictional heating was taken into account. This re- 
port discusses these pairs of solutions In greater detail and shows how 
the solutions are affected by systematic variations of the basic physi- 
cal parameters. It is also shown that a critical set of parametric 
values exists beyond which no fully developed solutions can be found. 


INTRODUCTION 

In recent analyses of combined forced- and natural- convection flow 
and heat transfer in enclosed regions. It was pointed out that the ef- 
fects of frictional heating could, in practical cases, be important 
(refs . 1 to 4) . Specifically, consideration was given to the fully de- 
veloped flow between two parallel surfaces which were subject to various 
thermal conditions and oriented parallel to the body force direction. 

As a result of retaining the frictional heating terms in the basic equa- 
tions, the final fourth-order ordinary differential equation was non- 
linear. The boundary value problems specified by the nonlinear equation 
and appropriate boundary conditions were at first solved approximately 
by an analytical iteration technique in references 1, 3, and 4. 

As a check of the accuracy of those solutions, the complete non- 
linear problem was solved by numerical integration on an IBM Card 
Programmed Calculator (CPC) . It was then discovered that the problems 
had either a pair of solutions or no solution, and examples of the pairs 
of solutions are given in references 1, 3, and 4. The solutions ob- 
tained analytically displayed no such behavior although they did approxi- 
mate one of each pair of solutions in limited ranges of the parameters. 
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A somewhat analogous situation was noted by Hartree (ref. 5) in 
his study of boundary-layer flows with adverse pressure gradients . The 
second solutions to that problem are presented and discussed further in 
references 6 and 7 . 

Since these results are unusual and interesting, and only a few 
examples are presented in the literature, the purpose of the present re- 
port is to discuss in some detail the pairs of solutions to the convec- 
tion problem. To this end, the method of obtaining these solutions is 
presented, and the influence of the various physical parameters on these 
solutions is discussed. 


BASIC EQUATIONS 

Consideration is given herein to the fully developed laminar flow 
of a viscous fluid between two vertical plane surfaces and subject to a 
vertical body force (see fig. l) . Fully developed flow means that the 
velocity components are independent of the axial distance. Therefore, 
the solutions obtained are valid away from the channel ends and, hence, 
apply to channels with large length- to-gap ratios. It is further as- 
sumed that the physical properties of the fluids are constants, except 
that the essential influence of density variations on the flow is ac- 
counted for by the introduction of the fluid volumetric expansion coef- 
ficient in the body force term. The latter procedure is given in detail 
in appendix B of reference 1 and is justified for liquids and for gases 
only if pressure differences are small relative to temperature differ- 
ences. The other influences of variable density and the variation of 
the expansion coefficient with temperature are also considered to be 
negligible. Fluids satisfying these assumptions will be called "quasi 
incompressible." 

Under the conditions stated it follows that there is only one non- 
zero velocity component and that it is a function only of the transverse 
coordinate Y. Further, the temperature can be expressed as the sum of 
a linear function of the vertical (axial) coordinate and an arbitrary 
function of the horizontal (transverse) coordinate as 

T* - AX + T(Y) (l) 

where A « dT * fbX is the axial temperature gradient. (All symbols are 
defined in appendix A.) The temperatures of the two surfaces are speci- 
fied either to be uniform or to vary linearly with the axial coordinate. 
In either case, the two surface temperatures can differ by a constant. 

The reduction of the basic continuity, momentum, and energy equa- 
tions by the above considerations is explicitly given in references 1, 

3, and 4. The equations expressing the conservation of momentum and en- 
ergy thus become 
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u" + T = CK (2) 

t" - (Ra)u + (u 1 ) 2 + aK = 0 ( 3 ) 

The continuity equation is satisfied identically. The terms (from left 
to right) in equation ( 2 ) denote the viscous, buoyancy, and axial pres- 
sure forces and in equation (3) represent the conduction, convection, 
frictional heating, and heat-source effects. The parameter C essen- 
tially specifies the temperature level or axial pressure gradient, and 
K is the frictional heating parameter. The symbol Ra denotes the 
Rayleigh number (product of Prandtl and Grashof numbers), and a is the 
internal-heat-source parameter. 


The Rayleigh n umb er is the natural- convection equivalent of the 
Peelet number (product of Prandtl and Reynolds numbers), and K is 
analogous to the Eckert number TJ^/ 

pation in forced flows. 


C P% 


which is a measure of the dissi- 


Elimination of T between equations (l) and ( 2 ) yields the afore- 
mentioned fourth- order nonlinear ordinary differential, equation: 


u iv + (Ra)u - (u 1 ) 2 - aK = 0 


( 4 ) 


Equation (4) applies for the three problems treated in references 1, 3, 
and 4. For the constant-wall- temperature case, Ra = 0 (see eq. (l))j 
heating from below is simulated by the case Ra < 0, whereas Ra > 0 is 
the conventional flow with linearly varying wall temperatures. 


BOUNDARY CONDITIONS 

The no-slip condition for viscous fluids is imposed in a 1 l cases, 
that is, 


u(0) = u(l) = 0 (5) 

In the first problem treated in reference 1, it was specified that 
the temperature is constant along each wall but that the two walls are 
not necessarily at the same temperature. Hence, A = 0 and the thermal 
boundary conditions for this case can, by use of equation ( 2 ), be 
written as 


u"(o) =■ (C - l)K 

(6a) 

u"(l) =» m(C - l)K 

(7a) 
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■where m, which is essentially a measure of the asymmetry of the thermal 
boundary conditions , is given by 

m - (c - e Vi /e V() )/(c - 1 ) (8a) 

The thermal boundary conditions for linearly varying wall tempera- 
tures are 


where 


u"(o) 

= CK 

(6b) 

u"(l) 

= mCK 

(7b) 

m = 1 - 

0w 1 /CAd 

(8b) 


The cases A > 0 and A < 0 (the latter corresponds to heating from 
below) are discussed in references 3 and 4, respectively. For the 
varying-wall- temperature cases the reference temperature in K and 
elsewhere (denoted by Sr for the purpose of unifying past work in the 
present report) is Ad and not & w . 


METHOD OF SOLUTION 

It can be seen from equations (2), (3), (6), and (7) that to specify 
a problem the parameters CK, Ra, aK, and m must be known. However , 
the solution of the problem is more easily obtained if the boundary con- 
dition on the second derivative at y = 1 is replaced by an arbitrary 
choice of the first derivative at y » 0. Then the third derivative at 
zero is varied until integration yields u(l) = 0. The resultant u"(l) 
yields the m to which the solution just obtained applies. The itera- 
tion method used in obtaining the solutions is discussed in detail in 
appendix B. 

The variation of m with u'(o) is represented on a curve called 
an m-curve. Such a curve is illustrated in the following schematic 
sketch: 
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In general, only the uppermost part of an m- curve is of practical inter- 
est, and only this part is presented in subsequent figures. 

Each point on an m-curve corresponds to a solution of the problem 
with the prescribed CK, Ra, and txK. The m-curves indicate that maxi- 
mum and minimum values of m exist for each CK, Ra, and crK beyond 
which no solutions of the type considered in the analysis (i.e., fully 
developed) exist and between which pairs of solutions exist. Thus, 
arbitrary values of the parameters cannot be specified a priori, and two 
distinct states of flow and heat transfer are indicated for each given 
set of conditions within the limits discussed above. Similar results 
are reported in reference 7 for boundary- layer flows. 

Because limits on parametric values were indicated, it became neces- 
sary to determine these limits in addition to computing velocity and tem- 
perature profiles . To this end, regions of the limit curves near the 
maximum m values were investigated in detail, and only those portions 
are included herein (i.e., see fig. 2). Some investigations were car- 
ried far enough to exhibit the general nature of these curves as dis- 
cussed previously. The minimum values of m indicated therefrom were 
generally very large negative numbers, and these were not considered to 
be of much interest from a practical viewpoint. 


RESUITS 

The dependence of the limit curves on the parameters K, Ra, and a 
can be discussed first. The constant C was arbitrarily taken to be 
zero for the constant-wall- temperature case (eqs. (6a) and (7a)) (essen- 
tially fixing the temperature level); for the varying-wall- temperature 
case, C was given the value -1 (eqs. (6b) and (7b)) so that the bound- 
ary conditions for all cases were the same. For a = 0 the value of 
the parameter CK is considered to be determined by a combination of 
individual C or K values. 

The effect of different K is shown in figure 2 plotted from table 
I for the ease Ra = a, = 0 (i.e., constant wall temperature and no in- 
ternal heat sources). The m-curves for other Ra and a, are of the 
same form. Recall that only parts of the m-curves are presented 
here. In general, the limiting value of m increases algebraically as 
K decreases algebraically. In other words, to increase the probability 
of a problem having a solution decrease the value of K in magnitude. 
Physically lower K values imply less frictional heating. On the 
basis of this observation and s imi lar behavior for Ra and a varia- 
tions (which also indicate the amount of heat input to the channel), it 
appears that the limit beyond which no solutions of the type considered 
herein exist depends on the heat input to the channel. The reason for 
this result is not as yet evident. However, another example where the 
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heat input is limited is in the familiar one- dimensional channel flow 
where thermal choking occurs. 

The effect of variation of the heat-source parameter a plotted 
from tables I and II is shown in figure 3. This shows that the greater 
the heat source, the less likely a solution exists for a particular 
problem. This trend is also in accord with the physical discussion 
given previously. 

The influence of Ra variation on the limit curves is shown in 
figure 4 plotted from table III. Increasing Ra makes it more likely 
that a particular problem has a solution. This, too, is in accord with 
the previous discussion, because larger Ra implies increased axial 
temperature gradients and, hence, a large transfer of heat along the 
channel axis . This axial heat transfer should permit the fluid to be 
subject to greater heat input. 

The negative Ra cases, simulating heating from below, were not 
studied in as great detail as those reported previously. However, re- 
sults similar to those for positive Ra were found for the cases In- 
vestigated and are discussed in some detail in reference 4. 

The m-eurves show that, if the effect of frictional heating is 
taken into account in internal combined forced- and natural-convection 
flows as considered herein, two unusual results occur which do not occur 
if frictional heating is neglected: (l) there exist critical sets of 

conditions beyond which no fully developed solutions exist and (2) where 
solutions exist, there are, in general, two solutions for each set of 
admissible parametric values. The first of these results has been dis- 
cussed, but the second one should be considered in more detail. Figure 
5 presents a set of representative pairs of velocity and temperature 
distributions taken from the previous work reported In reference 3 for 
Ra = 10, K = 10, a, = 0, and m = 2. These profiles are associated with 
the two circled points on the limit curve in figure 4. 

Since each point on a limit curve corresponds to a solution of a 
given problem, the effect of parametric variations on the flow velocities 
can be observed in figures 1 to 3. For example, larger velocities are 
obtained for increasing values of the frictional heating parameter K; 
also, for a given m, increasing the internal-heat-source parameter a 
or decreasing the Rayleigh number Ra causes one of the pair of flows 
to have greater velocities and the other lower velocities. Note the 
broadness of the m-curve for large Ra. This broadness indicates that 
the profiles of each velocity pair and of each temperature pair differ 
greatly In magnitude except for m near the maximum. 
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The flow with the superscript (l) agrees reasonably well with that 
computed by the analytic method previously mentioned. In general, the 
analytic solution agrees with the one of the pair that has the lower 
velocity extremes because the iteration procedure employed to obtain the 
analytic solution started with the no-frictional-heating case. Because 
less frictional heating is associated with the less intense solution, 
this solution is the one to be expected in practice. The other solution, 
however, denoted with the superscript (2), represents a flow with veloci- 
ties and temperature differences more than an order of magnitude greater 
than the first. The second solutions are connected intimately with the 
regenerative action of frictional heating in natural convection -wherein 
frictional heating acts as a heat source in the fluid and thereby leads 
to larger flow velocities. As an illustration of the role that the dis- 
sipation plays in these problems, consider the simplest case, namely, 
the one wherein the walls are kept at the same constant temperatures and 
no heat sources are in the fluid (cc = 0) . For this case, which is 
treated in reference 1, the less intense solution yields essentially the 
conduction (linear) temperature distribution. That is., no matter what 
kind of temperature profile exists in the fluid as it enters the channel, 
the profile becomes the conduction distribution by the time the flow has 
become fully developed because heat is extracted at the w alls in the 
entrance region (see sketch (a)). For the more intense solution 



the entering profile becomes parabolic and differs from the conduction 
profile considerably (see sketch (b) and fig. 5(b)). It is thus evi- 
dent that the frictional heating acts as a heat source in the fluid. 
Another significant fact that can be seen from figure 5 and the other 
examples of pairs of solutions given in references 1, 3, and 4 is that 
the mass flow and energy associated with the more intense solution are 



8 


NACA TN 4273 


greater than for the less intense solution. It, therefore, seems that, 
if the intense flow is to be obtained in practice, the flow would have 
to be started with the mass flow or energy at the higher level. ThiB 
might be possible, for example, by heating the fluid considerably as it 
enters the channel. 

The next question to arise in connection with the second solution 
concerns its stability; that is, since the frictional heat is dependent 
on the large velocities, would a small decrease in velocity, for example, 
decrease the frictional heating, and would this process keep repeating 
until the second solution attenuates. The stability of the symmetric 
case (m = l) has recently been analyzed, and the results, which will be 
published soon, show that this case is stable. 

Although pairs of solutions were obtained for three somewhat differ- 
ent problems, all the analyses were restricted to fully developed flows 
in which the fluid was only "quasi incompressible." Recently, S. Maslen 
studied the same problems as discussed herein, but he considered the 
fluid to be compressible and to have property variations. He treated 
the complete nonlinear problem by Galerkin's method, and not only did he 
determine explicitly the conditions under which the present analysis is 
valid, but Maslen too obtained pairs of solutions. These results are 
to be published soon. 


CONCLUDING REMARKS 

It can be stated from the work reported herein and from the work by 
Maslen that two flow and heat-transfer states are predicted theoretically 
for a fixed set of parametric values within certain limits for fully de- 
veloped combined forced- and natural- convection flow of certain real 
fluids. The only restriction of the analysis which has not been evaluated 
is that of assuming the flow to be fully developed. Definitive experi- 
ments will have to be performed to evaluate this limitation in order to 
establish definitely the physical existance of the second state of flow 
and heat transfer indicated herein. 

The analyses also show that a critical set of parameters exists 
beyond which no fully developed solutions can be found. Here, too, ex- 
perimental observation of what occurs for parametric values beyond the 
critical would be of importance. 

Lewis Flight Propulsion Laboratory 

National Advisory Committee for Aeronautics 
Cleveland, Ohio, March 20, 1958 
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SYMBOLS 

longitudinal (axial) temperature gradient 


constant, + p^f^ -y/(Pr ) d 4 /c p e R K 

specific heat at constant pressure 

characteristic length (specif ically, distance between plates) 
negative of X-component of body force per unit mass 
Grashof number, Bf^^d^/v^ 
step size 
(Pr)(Gr)Bf x d/c p 

thermal- conductivity coefficient 
constant defined by eq. (8) 
pressure 

Prandtl number, Cpjj/k 

heat due to heat sources 

Rayleigh number, (cpp/k) (pf-^Ad^/v 2 ) 

temperature function defined by eq. (l) 

temperature 

velocity 

reference velocity 
dimensionless velocity, U-y/(Pr)K/ 
perturbation function 
longitudinal (axial) coordinate 
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Y 

y 

a 


e 


transverse coordinate 
dimensionless transverse coordinate 
dimensionless heat-source parameter, Qd 2 /K0g 

volumetric expansion coefficient, 

temperature difference, T - T s for constant-wall- temperature 
case and T - T W q for linear -wall- temperature case 



0 R reference temperature, T W q - T s for constant-wall-temperature 

case and Ad for linear-wall-temperature case 

p absolute viscosity coefficient 

v kinematic viscosity coefficient 

p density 

p s for constant-wall- temperature case and p W Q for linear- 
wall- temperature case 

t dimensionless temperature differences, K (6/9%) 

Subscripts: 

s hydrostatic condition 

Wq wall conditions at y - 0 

w-j_ wall conditions at y » 1 

Superscript: 

Primes denote differentiation with respect to y 
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ITERATION TECHNIQUE 

Various integration methods were used: trapezoidal integration 

with end correction, three-point integration, five-point integration, 
Runge-Kutta integration, and finally Taylor’s series integration. All 
methods gave the same results to six significant figures, but the last 
named gave good approximate results with a step size of 0.2 and, hence, 
was used in most of the calculations. When detailed profiles were de- 
sired, integrations were repeated with a step size of 0.1 or 0.05. 

The Taylor series formulas were 
10-J 

u C<5)(y + h) =* u ( n+ j)y — for j =* 0, 1, 2, 3 

/ , n’ 

So 

where the superscript on u refers to the order of the derivative. 
Derivatives of u of higher order than the fourth were obtained by dif- 
ferentiating equation ( 4 ) . 

The usual steps in a straightforward iteration were: 

(1) Guess u"'(o). 

( 2 ) Integrate from. 0 to 1. 

( 3 ) Test whether u(l) is sufficiently close to zero. 

(a) If so, stop. 

(b) If not, choose a new u’"(o) and return to step ( 2 ). 

The only potentially bothersome step is (b) . After the first inte- 
gration, u”' ( 0 ) may be decreased if u(l) is positive and increased if 
u(l) is negative - but by how much is a matter of Judgment and experience. 
After two iterations, linear interpolation or extrapolation usually gives 
a good new u’" (0). Since an integration took 5 minutes, a machine oper- 
ator did not have great difficulty in keeping ahead of several cases run 
in rotation. However, in situations were only one m-curve was being 
studied, it was desirable to use a scheme which would guess a new u’”(o) 
automatically. Hence, the following perturbation problem was integrated 
along with the main problem: 
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Let u be the solution of the initial value problem: 

u iv “ (u') 2 " (Ra)u + aK 

u(0) « 0, u’(0) - b, u"(0) = CK, u ,M (0) a g (Bl) 

Let W be the solution of the initial. value problem: 

W iv = (w*) 2 - (Ra)W + ocK 

W(0) = 0, W* (0) « b, W"(0) = CK, W m (0) = g + e (B2) 

Let S = W - u. Then substituting in equation (B2) yields 
u iv + s iv _ u i + s«)2 _ Ra ( u + s) + ccK 

u(0) + S(o) > 0, u'(0) + S’(o) = b (B3) 

u"(0) + S"(0) =* CK, u m (0) + S ,ri (o) a d + e 
and simplifying by using equation (Bl) yields 

S iv = 2S'u' + (S') 2 - (Ra)S 

s(o) = 0, s’(0) a 0, S"(o) a 0, S m (0) a e (B4) 

Now let S = eZ, so that S' a eZ T , and so forth. Then substituting 
in equation (B4) and dividing by e yield 

Z iv = 2Z'u' + e(Z*) 2 - (Ra)Z 

Z(0) = Z'(0) = Z"(0) = 0, Z"'(0) => 1 (B5) 

As e approaches zero, the function Z approaches the function V, 
which is the solution of 

V iv a 2V'u' - (Ba)V 

V(0) - V’CO) => V"(0) a 0, V'" (0) = 1 (B6) 

From the preceding work, it can be observed that 

afrffTS) = lim - lirn ^ « lim Z(l) a V(l) (B7) 

e -> 0 e -*• 0 e -►O 
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Hence, a good approximation of a suggested change in g or u' 1 ' (o) 
is 

Ag= -u(l)/[du(l)/du«»(0)] = -u(l)/v(l) (B8) 

so that the formula after any trial integration for the new u 1 " (o) is 
new u M, (0) = old u n, (o) - u(l)/v(l). 

Because of storage limitations of 38 ten-digit numbers, only four 
terms of the Taylor series for V were used, but even this yielded good 
convergence to a result in several trials . General-purpose floating 
point control panels were used for all calculations. 
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Figure 1. - Schematic sketch of configuration considered. 
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Figure 3. - The m-curves for various values of heat- 
source parameter a. Ra = 0; K = 10. 
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(a) Velocity. 


Figure 5. - A representative pa_Ir of profiles for Ra = 10, 
K = 10, m * 2, a = 0, and C - -1. 







